Imports

In [1]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from mpl_toolkits.mplot3d import Axes3D
import plotly
import plotly.express as px
from matplotlib import animation
from IPython.display import clear_output
from life_saving_tools.maths.ode import ODE
from IPython.display import Video, HTML
plt.rcdefaults()
o = ODE()
plotly.offline.init_notebook_mode()

Single-Particle Motions

Making Some Functions

Here, I'll define some functions that will be used in the following sections.

In [2]:
def solve_ode(fns, ts, ics, modify=False):
    """
    A function to solve an ODE using the `ODE` class's `solve` method.

    Parameters
    ----------
    fns : list
        A list of functions that represent the ODE's differential equations.
    ts : list
        A list of times at which the ODE will be solved.
    ics : tuple
        A tuple of initial conditions for the ODE.
    modify : bool, optional
        A boolean that determines whether the ODE will be solved with the
        modification of the initial conditions. The default is `False`.

    Returns
    -------
    np.array:
        A numpy array of the solutions to the ODE.
    """
    t0, x0 = ics
    o = ODE()
    xs = np.zeros((len(ts), len(fns)))
    for i in range(len(ts)):
        if modify:
            if i<=10:
                xs[i] = o.solve(funcs = fns, x0 = t0, y0=x0, x=ts[i])
            else:
                xs[i] = o.solve(funcs = fns, x0 = np.array([ts[i-5]]*len(fns)), y0=xs[i-5], x=ts[i])
        else:
            xs[i] = o.solve(funcs = fns, x0 = t0, y0=x0, x=ts[i])
    return xs
In [3]:
def show_results(xs, ts, title="", xlabel="Time", ylabel="Values", **kwargs):
    """
    Plots all the columns of xs with time

    Parameters
    ----------
    xs : np.array
        A numpy array of the solutions to the ODE.
    ts : np.array
        A numpy array of the times at which the ODE was solved.
    title : str
        The title of the plot.
    xlabel : str
        The label for the x-axis.
    ylabel : str
        The label for the y-axis.
    **kwargs : dict
        A dictionary of keyword arguments to be passed to the `plt.plot`.

    Returns
    -------
    None
    """
    plt.figure(figsize=(10,6))
    for i in range(len(xs[0, :])):
        plt.plot(ts, xs[:, i], label=f"$u_{i+1}$", **kwargs)
    plt.xlim(0, ts[-1])
    plt.legend()
    plt.title(title)
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.show()
In [4]:
def plot_x_with_y(xs, t=None, title="", **kwargs):
    """
    Plots x with y using plotly.

    Parameters
    ----------
    xs : np.array
        A numpy array of the solutions to the ODE.
    title : str
        The title of the plot.
    **kwargs : dict
        A dictionary of keyword arguments to be passed to the `plotly.express.line`.

    Returns
    -------
    None
    """
    color = t if t is not None else np.random.rand(3)
    fig = px.scatter(
        x=xs[:, 0],
        y=xs[:, 1],
        color=color,
    )
    fig.update_layout(
        title=title,
        xaxis_title="x",
        yaxis_title="y",
        width=600,
        height=600,
        **kwargs
    )
    fig.show()
In [5]:
def plot_x_with_y_and_z(xs, zs, title="", **kwargs):
    """
    Plots x with y and z using plotly.

    Parameters
    ----------
    xs : np.array
        A numpy array of the solutions to the ODE.
    zs : np.array
        A numpy array of the solutions to the ODE containg the z-coordinate.
    title : str
        The title of the plot.
    **kwargs : dict
        A dictionary of keyword arguments to be passed to the `plotly.express.line`.
    
    Returns
    -------
    None
    """
    fig = px.line_3d(x=xs[:, 0], y=xs[:, 1], z=zs)
    fig.update_layout(
        title=title,
        xaxis_title="x",
        yaxis_title="y",
        width=600,
        height=600,
        **kwargs
    )
    fig.show()
In [6]:
def animate_trajectory(
    xs, t, title="", num_points_per_time_period=100, interval=20, fps=30, **kwargs
):
    """
    Animates the trajectory of the particle using matplotlib.

    Parameters
    ----------
    xs : np.array
        A numpy array of the solutions to the ODE.
    t : np.array
        A numpy array of the times at which the ODE was solved.
    title : str
        The title of the plot.
    num_points_per_time_period : int
        The number of points to be plotted per time period. default is 100.
    **kwargs : dict
        A dictionary of keyword arguments to be passed to the `plt.figure`.
    """
    fig = plt.figure(figsize=(8, 8), **kwargs)
    # Setting axes limits
    margin = 0.2
    xlims = [xs[:, 0].min() - margin, xs[:, 0].max() + margin]
    ylims = [xs[:, 1].min() - margin, xs[:, 1].max() + margin]
    axis = plt.axes(xlim=xlims, ylim=ylims)
    # Setting labels for the axes
    plt.xlabel("x")
    plt.ylabel("y")
    plt.title(title)

    (line,) = axis.plot([], [], lw=2)

    def get_position(i):
        return xs[i, 0], xs[i, 1]

    # initializing the x and y data holders
    def init():
        line.set_data([], [])
        return (line,)

    # initializing empty values for x and y
    xdata, ydata = [], []
    legend = plt.legend()

    # animation function
    def animate(i):
        if (i + 1) % 50 == 0:
            clear_output(wait=True)
            print(f"TimeStamp: {i+1} of {len(t)}.")
        x, y = get_position(i)
        # appending values to the previously initialized x and y data holders
        xdata.append(x)
        ydata.append(y)
        line.set_label("t = " + str(round(t[i], 2)))
        # Changing the line color
        colors = ["r", "g", "b", "m", "k"]
        j = i // num_points_per_time_period
        j = j % len(colors)
        line.set_color(colors[j])

        legend = plt.legend(loc="upper right")
        line.set_data(xdata, ydata)

        return line, legend

    # calling the animation function
    anim = animation.FuncAnimation(
        fig, animate, init_func=init, frames=len(t), interval=interval, blit=True
    )

    # saves the animation
    file_name = title.replace(" ", "_") + ".mp4"
    anim.save(file_name, writer="ffmpeg", fps=fps)
    return file_name

Uniform B and Uniform (or time varying) E

E=0

Solving Numerically

Define the following variables: $$ \begin{align*} x &= u_1\\ y &= u_2\\ \dot{x} &= u_3\\ \dot{y} &= u_4\\ \end{align*} $$ Using this notation, we can write the equation of motion for the particle: $$ \begin{align*} \dot{u_1} &= u_3\\ \dot{u_2} &= u_4\\ \dot{u_3} &= c u_4\\ \dot{u_4} &= -c u_3\\ \end{align*} $$ Here c is a constant which is equal to $\frac{qB}{m}$.

The z component of the particle's velocity is a constant. That is, for z component, we have: $$ z = v_zt+z_0 $$

We can solve the system of above equation as following.

In [7]:
#Defining the constants
c = 1
vz = 1
z0=0

#Defining the function to be solved
f1 = lambda t, u1, u2, u3, u4: u3
f2 = lambda t, u1, u2, u3, u4: u4
f3 = lambda t, u1, u2, u3, u4: c*u4
f4 = lambda t, u1, u2, u3, u4: -c*u3
funcs = [f1, f2, f3, f4]

#The time period
T = 2*np.pi/c
# Number of time periods to plot
NUM_TIME = 1
num_points_per_time_period = 100
T_last = NUM_TIME*T + T*0.01
t = np.linspace(0, T_last, NUM_TIME*num_points_per_time_period) #Plotting till `NUM_TIME` period

# Initial conditions
t0 = np.array([0, 0, 0, 0])
x0 = np.array([0, 0, 1, 1])

# Solving the ODE for x and y
xs = solve_ode(funcs, t, (t0, x0), modify=False)
# Getting z
zs = vz*t + z0

# Plotting the Results
title="Partical in Constant B and zero E"
show_results(xs, t, title)

Here is a plot showing the motion of the particle in the x-y plane.

In [8]:
#Plotting x with y
plot_x_with_y(xs,t, title)

And here is a 3d plot.

In [9]:
#Plotting x with y and z
plot_x_with_y_and_z(xs, zs, title)

Using the Analytical Solution

The analytical solution is: $$ \begin{align*} x &= x_0+r_l\cos(\omega_c t)\\ y &= y_0+r_l\sin(\omega_c t)\\ z &= z_0+vt \end{align*} $$ where: $$ \begin{align*} r_l &= \frac{mv}{qB}\\ \omega_c &= \frac{qB}{m} \end{align*} $$

In [10]:
x0 = 1
y0 = 1
ts = np.linspace(0, 10, 100)
m = 1
v=1
q=1
B=1
rl = (m*v)/(q*B)
omega = (q*B/m)

xs_a = x0 + rl*np.sin(omega*ts)
ys_a = y0 + rl*np.cos(omega*ts)

plt.figure(figsize=(8, 8))
plt.plot(xs_a, ys_a, 'b-')
plt.scatter(x0, y0, c='r', s=100, label="Center")
plt.xlabel('x')
plt.ylabel('y')
plt.title('Single-Particle Motion')
plt.legend()
plt.show()

Animate the Trajectory

An animation for one time period.

In [ ]:
file = animate_trajectory(xs, t, title=title, fps=10)
In [12]:
Video(file)
Out[12]:

5 Time Periods

In [ ]:
#The time period
T = 2*np.pi/c
# Number of time periods to plot
NUM_TIME = 5
num_points_per_time_period = 100
T_last = NUM_TIME*T + T*0.01
t = np.linspace(0, T_last, NUM_TIME*num_points_per_time_period) #Plotting till `NUM_TIME` period

# Initial conditions
t0 = np.array([0, 0, 0, 0])
x0 = np.array([0, 0, 1, 1])

# Solving the ODE for x and y
xs = solve_ode(funcs, t, (t0, x0), modify=False)
title = "Partical in Constant B and zero E for 5 T"
file = animate_trajectory(xs, t, title=title)
In [14]:
Video(file)
Out[14]:

Simulating for 5 time periods with modified initial conditions.

In [ ]:
#The time period
T = 2*np.pi/c
# Number of time periods to plot
NUM_TIME = 5
num_points_per_time_period = 100
T_last = NUM_TIME*T + T*0.01
t = np.linspace(0, T_last, NUM_TIME*num_points_per_time_period) #Plotting till `NUM_TIME` period

# Initial conditions
t0 = np.array([0, 0, 0, 0])
x0 = np.array([0, 0, 1, 1])

# Solving the ODE for x and y
xs = solve_ode(funcs, t, (t0, x0), modify=True)
title = "Partical in Constant B and zero E for 5 T modified"
file = animate_trajectory(xs, t, title=title)
In [16]:
Video(file)
Out[16]:

Finite Constant E

Solving Numerically

In case of finite E, we get the following equation of motion: $$ \begin{align*} \dot{v_x} &= \frac{qE_x}{m}+w_c v_y\\ \dot{v_y} &= -w_c v_x\\ v_z &= \frac{qE_z}{m}t + v_{z_0} \end{align*} $$

The solution to the transverse part can be determined by defining the following variables: $$ \begin{align*} x &= u_1\\ y &= u_2\\ \dot{x} &= u_3\\ \dot{y} &= u_4\\ \end{align*} $$ Using this, the equation of motion for the particle in the x-y plane is: $$ \begin{align*} \dot{u_1} &= u_3\\ \dot{u_2} &= u_4\\ \dot{u_3} &= \frac{qE_x}{m}+ \omega_c u_4\\ \dot{u_4} &= -\omega_c u_3\\ \end{align*} $$ Let's solve this.

In [17]:
def Ex(t):
    #Using a constant value
    return 1

def Ez(t):
    #Using a constant value
    return 1
In [18]:
q, m, omega = 1, 1, 1
f1 = lambda t, u1, u2, u3, u4: u3
f2 = lambda t, u1, u2, u3, u4: u4
f3 = lambda t, u1, u2, u3, u4: q * Ex(t) / m + omega * u4
f4 = lambda t, u1, u2, u3, u4: -omega * u3
funcs = [f1, f2, f3, f4]
# The time period
T = 2 * np.pi / omega
# Number of time periods to plot
NUM_TIME = 5
T_last = NUM_TIME * T + T * 0.01
t = np.linspace(0, T_last, NUM_TIME * num_points_per_time_period)  # Plotting till one time period

# Initial conditions
t0 = np.array([0, 0, 0, 0])
x0 = np.array([0, 0, 1, 1])

# Solving the ODE for x and y
xs = solve_ode(funcs, t, (t0, x0), modify=False)

# Plotting the Results
title="Partical in Constant B and constant E"
show_results(xs, t, title)
In [19]:
plot_x_with_y(xs, t,title)

To make a 3d plot, let's solve for the z-coordinate too.

The equation of motion for the particle in the z axis is:$$v_z = \frac{qE_z}{m}t + v_{z_0}$$

To solve this, we define: $$ \begin{align*} z &= u_1\\ \dot{z} &= u_2\\ \end{align*} $$ Using this, the equation of motion becomes: $$ \begin{align*} \dot{u_1} &= u_2\\ \dot{u_2} &= \frac{qE_z}{m}t + v_{z_0} \end{align*}

Let's solve this:

In [20]:
vz0 = 0
f1z = lambda t, u1, u2: u2
f2z = lambda t, u1, u2: q*Ez(t)*t/m + vz0
funcsz = [f1z, f2z]
#The time period
T = 2*np.pi/omega
# Number of time periods to plot
NUM_TIME = 5
T_last = NUM_TIME*T + T*0.01
t = np.linspace(0, T_last, NUM_TIME*num_points_per_time_period) #Plotting till NUM_TIME time period
t0 = np.array([0, 0])
x0 = np.array([0, 1])
zs = solve_ode(funcsz, t, (t0, x0), modify=False)
title = "Partical in Constant B and constant E"
show_results(zs, t, title)
In [21]:
plot_x_with_y_and_z(xs, zs[:, 0], title)

Animate the Trajectory

In [ ]:
title = "Partical in Constant B and constant E for 5 T"
file = animate_trajectory(xs, t, title=title)
In [23]:
Video(file)
Out[23]:
In [ ]:
NUM_TIME = 5
T_last = NUM_TIME * T + T * 0.01
t = np.linspace(0, T_last, NUM_TIME * num_points_per_time_period)  # Plotting till one time period

# Initial conditions
t0 = np.array([0, 0, 0, 0])
x0 = np.array([0, 0, 1, 1])

# Solving the ODE for x and y
xs = solve_ode(funcs, t, (t0, x0), modify=True)

title = "Partical in Constant B and constant E for 5 T modified"
file = animate_trajectory(xs, t, title=title)
In [25]:
Video(file)
Out[25]:

Let's reapeat the same simulation but with a varying E.

Finite Time Varying E

We'll assume that the electric filed is varying sinusoidally with time in x direction and is zero otherwise. The equation of motion remains the same.

In [26]:
q, m, omega = 1, 1, 1
def Ex(t, omega=omega):
    #Using a constant value
    return np.sin(omega * t)

def Ez(t):
    #Using a constant value
    return 0
In [ ]:
NUM_TIME = 5
T_last = NUM_TIME * T + T * 0.01
t = np.linspace(0, T_last, NUM_TIME * num_points_per_time_period)  # Plotting till one time period

# Initial conditions
t0 = np.array([0, 0, 0, 0])
x0 = np.array([0, 0, 1, 1])

# Solving the ODE for x and y
xs = solve_ode(funcs, t, (t0, x0), modify=False)

title = "Partical in Constant B and time varying E for 5 T"
file = animate_trajectory(xs, t, title=title)
In [28]:
plot_x_with_y(xs, t, title)
In [29]:
Video(file)
Out[29]:
In [30]:
zs = vz*t + z0
plot_x_with_y_and_z(xs, zs, title)
In [31]:
xs_modified = solve_ode(funcs, t, (t0, x0), modify=True)
plot_x_with_y_and_z(xs_modified, zs, title=title+" modified")
In [ ]:
file = animate_trajectory(xs_modified, t, title=title)
In [33]:
Video(file)
Out[33]:

Let's reapeat the same simulation but with a E varying wit x.

Finite Spatial Varying E

We'll assume that the electric filed is varying sinusoidally with time in x direction and is zero otherwise. The equation of motion remains the same.

In [34]:
q, m, omega = 1, 1, 1
def Ex(t, omega=omega):
    #Using a constant value
    return np.sin(omega * t)

def Ez(t):
    #Using a constant value
    return 0
In [35]:
q, m, omega = 1, 1, 1
f1 = lambda t, u1, u2, u3, u4: u3
f2 = lambda t, u1, u2, u3, u4: u4
f3 = lambda t, u1, u2, u3, u4: q * Ex(u1) / m + omega * u4
f4 = lambda t, u1, u2, u3, u4: -omega * u3
funcs = [f1, f2, f3, f4]
# The time period
T = 2 * np.pi / omega
NUM_TIME = 5
T_last = NUM_TIME * T + T * 0.01
t = np.linspace(0, T_last, NUM_TIME * num_points_per_time_period)  # Plotting till one time period

# Initial conditions
t0 = np.array([0, 0, 0, 0])
x0 = np.array([0, 0, 1, 1])

# Solving the ODE for x and y
xs = solve_ode(funcs, t, (t0, x0), modify=True)

# Plotting the Results
title="Partical in Constant B and spatial varying E"
show_results(xs, t, title)
In [36]:
plot_x_with_y(xs, t, title)

The transverse component of the equation of motion also remains the same.

In [37]:
vz0 = 0
f1z = lambda t, u1, u2: u2
f2z = lambda t, u1, u2: q*Ez(t)*t/m + vz0
funcsz = [f1z, f2z]
#The time period
T = 2*np.pi/omega
# Number of time periods to plot
NUM_TIME = 5
T_last = NUM_TIME*T + T*0.01
t = np.linspace(0, T_last, NUM_TIME*num_points_per_time_period) #Plotting till NUM_TIME time period
t0 = np.array([0, 0])
x0 = np.array([0, 1])
zs = solve_ode(funcsz, t, (t0, x0), modify=False)
plot_x_with_y_and_z(xs, zs[:, 0], title)
In [ ]:
file = animate_trajectory(xs, t, title=title)
In [39]:
Video(file)
Out[39]:

A Generalized Implementation

The Generalization

Here, we'll implement the equations of motion for a particle in a general electric and magnetic field. The equations of motion are: $$ \begin{align*} m\dot{v_x} &= v_yB_z-v_zB_y+qE_x\\ m\dot{v_y} &= v_zB_x-v_xB_z+qE_y\\ m\dot{v_z} &= v_xB_y-v_yB_x+qE_z\\ \end{align*} $$

Same as before, we can define the variables: $$ \begin{align*} x &= u_1\\ y &= u_2\\ z &= u_3\\ \dot{x} &= u_4\\ \dot{y} &= u_5\\ \dot{z} &= u_6\\ \end{align*} $$

Using these, the equations of motion are: $$ \begin{align*} \dot{u_1} &= u_4\\ \dot{u_2} &= u_5\\ \dot{u_3} &= u_6\\ \dot{u_4} &= \frac{qB_z}{m}u_5- \frac{B_yu_6}{m}+\frac{qE_x}{m}\\ \dot{u_5} &= \frac{qB_x}{m}u_6- \frac{B_zu_4}{m}+\frac{qE_y}{m}\\ \dot{u_6} &= \frac{qB_y}{m}u_4- \frac{B_xu_5}{m}+\frac{qE_z}{m} \end{align*} $$

Let's solve this.

In [48]:
#Defining E and B
Ex = lambda x, y, z, t: 0
Ey = lambda x, y, z, t: 0
Ez = lambda x, y, z, t: 0

Bx = lambda x, y, z, t: 0
By = lambda x, y, z, t: 0
Bz = lambda x, y, z, t: 1

#Defining the constants
q, m = 1, 1
In [49]:
#Defining the functions for the ODE
u1 = lambda t, u1, u2, u3, u4, u5, u6: u4
u2 = lambda t, u1, u2, u3, u4, u5, u6: u5
u3 = lambda t, u1, u2, u3, u4, u5, u6: u6
u4 = lambda t, u1, u2, u3, u4, u5, u6: (q * Bz(u1, u2, u3, t) / m) * u5 - (q * By(u1, u2, u3, t) / m) * u6 + (q * Ex(u1, u2, u3, t) / m)
u5 = lambda t, u1, u2, u3, u4, u5, u6: (q * Bx(u1, u2, u3, t) / m) * u6 - (q * Bz(u1, u2, u3, t) / m) * u4 + (q * Ey(u1, u2, u3, t) / m)
u6 = lambda t, u1, u2, u3, u4, u5, u6: (q * By(u1, u2, u3, t) / m) * u4 - (q * Bx(u1, u2, u3, t) / m) * u5 + (q * Ez(u1, u2, u3, t) / m)

functions = [u1, u2, u3, u4, u5, u6]
In [50]:
# The time period
omega = 1
T = 2 * np.pi / omega
# Number of time periods to plot
NUM_TIME = 5
T_last = NUM_TIME * T + T * 0.01
t = np.linspace(0, T_last, NUM_TIME * num_points_per_time_period)  # Plotting till NUM_TIME period

# Initial conditions
t0 = np.array([0, 0, 0, 0, 0, 0])
x0 = np.array([0, 0, 0, 1, 1, 1])

# Solving the ODE for x and y
xs = solve_ode(functions, t, (t0, x0), modify=True)

show_results(xs, t, "Particle in a general field.")
In [52]:
plot_x_with_y(xs, t, "Particle in a general field")
In [ ]:
plot_x_with_y_and_z(xs, xs[:, 2], "Particle in a general field")
In [88]:
t0 = np.array([0, 0, 0, 0, 0, 0])
x0 = np.array([0, 0, 0, 1, 1, 1])
ic = (t0, x0)
def general_fields(Es, Bs, time_periods=5, ic = ic, title="Particle in a general field", position=False, modify=True):
    Ex, Ey, Ez = Es
    Bx, By, Bz = Bs
    #Defining the constants
    q, m = 1, 1
    #Defining the functions for the ODE
    u1 = lambda t, u1, u2, u3, u4, u5, u6: u4
    u2 = lambda t, u1, u2, u3, u4, u5, u6: u5
    u3 = lambda t, u1, u2, u3, u4, u5, u6: u6
    u4 = lambda t, u1, u2, u3, u4, u5, u6: (q * Bz(u1, u2, u3, t) / m) * u5 - (q * By(u1, u2, u3, t) / m) * u6 + (q * Ex(u1, u2, u3, t) / m)
    u5 = lambda t, u1, u2, u3, u4, u5, u6: (q * Bx(u1, u2, u3, t) / m) * u6 - (q * Bz(u1, u2, u3, t) / m) * u4 + (q * Ey(u1, u2, u3, t) / m)
    u6 = lambda t, u1, u2, u3, u4, u5, u6: (q * By(u1, u2, u3, t) / m) * u4 - (q * Bx(u1, u2, u3, t) / m) * u5 + (q * Ez(u1, u2, u3, t) / m)

    functions = [u1, u2, u3, u4, u5, u6]

    # The time period
    omega = 1
    T = 2 * np.pi / omega
    # Number of time periods to plot
    NUM_TIME = time_periods
    T_last = NUM_TIME * T + T * 0.01
    t = np.linspace(0, T_last, NUM_TIME * num_points_per_time_period)  # Plotting till NUM_TIME period

    # Initial conditions
    t0, x0 = ic

    # Solving the ODE for x and y
    xs = solve_ode(functions, t, (t0, x0), modify=modify)

    #plot
    plot_x_with_y(xs, t, title)
    plot_x_with_y_and_z(xs, xs[:, 2], title)
    if position:
        show_results(xs, t, title)
    return xs, t

Using the Generalized Equations

Constant B and E = 0

In [114]:
#Defining E and B
Ex = lambda x, y, z, t: 0
Ey = lambda x, y, z, t: 0
Ez = lambda x, y, z, t: 0

Bx = lambda x, y, z, t: 0
By = lambda x, y, z, t: 0
Bz = lambda x, y, z, t: 1

Es = [Ex, Ey, Ez]
Bs = [Bx, By, Bz]
xs, t = general_fields(Es, Bs, 5, title = "Particle in a general field", position=True)

B = 0 and E depending on position

In [121]:
#Defining E and B
Ex = lambda x, y, z, t: 0
Ey = lambda x, y, z, t: -y
Ez = lambda x, y, z, t: 0

Bx = lambda x, y, z, t: 0
By = lambda x, y, z, t: 0
Bz = lambda x, y, z, t: 0

Es = [Ex, Ey, Ez]
Bs = [Bx, By, Bz]
t0 = np.array([0, 0, 0, 0, 0, 0])
x0 = np.array([0, 0, 0, 1, 1, 0])
ic = (t0, x0)
xs, t = general_fields(Es, Bs, 5,ic, "Particle in a general field", position=True)
In [124]:
px.scatter(x = xs[:, 4], y = xs[:, 1], height=600, width=600)